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Abstract 

The method of space-time conservation element and solution element is a nontra- 
ditional numerical method designed from a physicist’s perspective, i.e., its development 
is based more on physics than numerics. It uses only the simplest approximation tech- 
niques and yet is capable of generating nearly perfect solutions for a 2-D shock reflection 
problem used by Helen Yee and others. In addition to providing an overall view of the 
new method, in this paper we shall introduce a new concept in the design of implicit 
schemes, and use it to construct a highly accurate solver for a convection-diffusion 
equation. It will be shown that, in the inviscid case, this new scheme becomes explicit 
and its amplification factots are identical to those of the Leapfrog scheme. On the 
other hand, in the pure diffusion case, its principal amplification factor becomes the 
amplification factor of the Crank-Nicolaon scheme . 



Figure 1 . — Pressure distribution for a shock 
reflection problem (a flow of Mach number 2.9 
enters from the left; the shock is reflected by a 
wall on the back). 

1. Introduction 

The method of space-time conservation element and solution element is a new 

numerical discretization method for solving conservation laws [1-11]. It is designed to 



overcome several key limitations of the well-established methods-i.e., finite difference, 
finite volume, finite element, and spectral methods. At the root of its development 
is a constant drive toward simplicity, generality, and accuracy. To appreciate the 
significance of this effort, consider a 2-D time-marching Euler solver developed using the 
new method [2,10]. Even though it does not use (i) any approximation techniques that 
are more complicated than Taylor’s expansion, (ii) any mesh-refinement techniques, 
(iii) any monotonicity constraints, (iv) any characteristics-based techniques, or (v) 
any ad hoc techniques that are used only in the neighborhood of a discontinuity, this 
solver is capable of generating highly accurate solutions for a 2-D shock reflection 
problem used by Helen Yee and others [12]. As shown in Fig. 1, both the incident 
and the reflected shocks can be resolved by a single data pant without the presence of 
numerical oscillations near the discontinuity. 

The Introduction section of [1] begins with a lengthy discussion that focuses more 
on physics than on numerics. From this discussion, readers can understand (i) the 
considerations that motivate the new method, and (ii) the key differences that separate 
it from the traditional methods mentioned above. 

By using a set of design principles that are extracted from the above discussion, 
several two-level explicit schemes were constructed in [1,8] to solve (i) the pure convec- 
tion equation 

dujdt + odu/dz = 0 (1-1) 

and (ii) the convection-diffusion equation 

dufdt + adu/dx — p&ufdx 1 = 0 (1.2) 

where the convection velocity a, and the viscosity coefficient p (> 0) are constants. 
These schemes were then extended to solve the 1-D time-dependent Euler and Navier- 
Stokes equations of a perfect gas [1,8]. Moreover, except for the Navier-Stokes solver, 
the above 1-D schemes have been generalized to their 2-D counterparts [2,10]. Because 
of the inherent simplicity and generality of the current method, the above multidimen- 
sional generalization is a straightforward matter. Also, as a result of the similarity in 
their designs, each of the above 2-D schemes shares with its 1-D version virtually the 
same fundamental characteristics. 

In addition to providing an overall view of the new method, in this paper we shall 
describe a simple and innovative approach by which accurate implicit time-marching 
solvers can be constructed using the new method. A striking feature of this new 
treatment is that the modeling of the diffusion-related terms involves interpolation 
between neighboring mesh points while that of the convection-related term does not. 
As a preliminary, first we shall discuss the pros and cons of explicit and implicit schemes. 



Figure 2. — (a) An initial-value problem, (b) An initial-value/boundary-value problem. 


For a two-level explicit scheme, the value of a solution at any mesh point has a 
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finite domain of dependence at the previous time level. As an example, consider a 
finite-difference solver for Eq. (1.1). Let ti”, the mesh value of ti at any mesh point 
(j,n) (point P in Fig. 2(a)) be determined by u?Zi , ti" -1 , and «"+*. Then the domain 
of dependence of uj at the (n — l)th time level contains three mesh points. Also one 
can see that u" is dependent only on the initial data given on the line segment AB. 

For an initial-value problem, such as a time-dependent Euler problem or a problem 
involving Eq. (1.1), the solution at any point in space-time also has a finite domain of 
dependence on the initial plane. As a result, explicit schemes could be ideal solvers for 
such a problem if they satisfy the requirement that the physical domain of dependence 
be a subset of the numerical domain of dependence. 

On the other hand, the solution of an initial-value/boundary-value problem at 
any point in space-time is dependent on the initial data and the boundary data up 
to the timf of the point under consideration. As an example, consider a problem 
involving Eq. (1.2). As shown in Fig. 2(b), the solution at point P is dependent on the 
initial/boundary data given on A' B, BC , and CD 1 where A! , P , and D are at the same 
tim» level. Let this problem be solved using the explicit scheme that was explained 
lining Fig. 2(a). Let P also be a mesh point O', n). Then uj is dependent only on the 
initial/boundary data given on AB, BC and CD. It is completely independent of those 
data given on AA' and DD. Contrarily, if the same problem is solved using an implicit 
scheme, then ti” is dependent on the initial/boundary data given on A'B, BC, and 
CD 1 . In other words, the numerical domain of dependence of the implicit scheme is 
consistent with the physical domain of dependence of the problem under consideration. 
Two observations can be made as a result of the above discussions. 

(a) Generally an explicit scheme is not an ideal solver for an initial-value/boundary- 
value problem. Because a time-dependent Navier-Stokes problem is such a prob- 
lem, the above argument implies that an explicit scheme cannot be used to solve 
a time-dependent Navier-Stokes problem except for the special circumstance in 
which errors caused by neglecting certain initial/boundary data (such as those 
given on AA' and DD* in Fig. 2(a)) are relatively small. The factors that help 
achieve the above special circumstance include: (i) a small time-step size to spatial- 
mesh interval ratio, (ii) a small time rate of change of boundary data, and (iii) 
a small contribution of the viscous terms in the Navier-Stokes equations relative 
to that of the inertial terms. Note that condition (iii) may be met by a high- 
Reynolds- number flow. 

(b) Generally an implicit scheme is not an ideal solver for an initial-value problem. 
This is because the domain of dependence of the former is far greater than the 
latter and, as a result, an implicit solution tends to be contaminated by extraneous 
information. 

This completes the discussion of explicit and implicit schemes. A review of an 
explicit scheme [1,8] will precede the construction of an implicit solver. 

2. Review of an Explicit Scheme 

Let Eq. (1.1) be in a dimensionless form. Let x\ — x and xi = t be considered as 
the coordinates of a two-dimensional Euclidean space f? 2 - By using Gauss’ divergence 
theorem in the space-time f? 2 , it can be shown that Eq. (1.1) is the differential form of 
the integral conservation law 

<f h • ds = 0 (2.1) 

/s(V) 

Here (i) S(V) is the boundary of an arbitrary space-time region V in Ei, (ii) h = (au, ti) 
is a current density vector in Ei, and (iii) ds — don with do and n, respectively, being 
the area and the outward unit normal of a surface element on S(V). Note that (i) 
h-dsis the space-time flux of h leaving the region V through the surface element ds, 
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and (ii) all mathematical operations can be carried out as though f?a were an ordinary 
two-dimensional Euclidean space. 

At this juncture, note that the conservation law Eq. (2.1) appears in a form in 
which space and time are unified and treated cm the same footing. This unity of space 
and tim e is a key characteristic that distinguishes the current method from most of the 
traditional methods. 
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Figure 3. — The SEs and CEs. (a) The relative positions 
of SEs and CEs. (b) SEQ, n). (c) CE_Q, n). (d) CE+Q, n). 


In the following, we shall review the inviscid version of the explicit a-p scheme 
[1,8]. To achieve consistency of the notations used in this and the next sections, the 
notations used here will be slightly different from those used in [1,8]. 

Let fli denote the set of mesh points (j, n) in E? (dots in Fig. 3(a)) where n = 
0,±1, ±2, ±3, . . ., and, for each n, j=n± 1, n±3,n±5, . . .. There is a solution element 
(SE) associated with each (j, n) € fii- Let the solution element SE (j,n) be the space- 
time region bounded by the dashed curve depicted in Fig. 3(b). It includes a horizontal 
line segment, a vertical line segment, and their immediate neighborhood. 

For any (x,t) € SE(j, n), u(z,t), and h(x,t), respectively, are approximated by 
t i*(x,f ; j,n) and h*(x,t ; j,n) which we shall define shortly. Let 

ti* (x, t ; j, n) = u? + (u,)j?(x - Xj ) + (ti, )”(* -t n ) (2.2) 


where (i) ti”, (ti*)”, and (a*)" are constants in SE(j, n), and (ii) (xj,t") are the coor- 
dinates of the mesh point ( j , n). Note that 


*■(*.*■*»)-« =M7 (M ) 


Moreover, if we identify u", (ti*)”, and (u«)", respectively, with the values of u, du/dx, 
and du/dt at (xj,t n ), the expression on the right side of Eq. (2.2) becomes the first- 
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order Taylor’s expansion of u(x,t) at (x Thus, u", («*)", and (u «)" are the 
numerical analogues of the values of u, du/dx, and du/dt at (x,-,*"), respectively. 

We shall require that u = u*(x, t ; j, n) satisfy Eq. (1.1) within SE(j, n), i.e., 

(u ,)? = -«(«,)? (2-4) 

Combining Eqs. (2.2) and (2.4), one has 

u*(x, t ,j, n) = u" 4- («r)" [(* - x,) - a (f - < n )J (x,<) 6 SE(j, n) (2.5) 


Because h — (au, u), we define 

h*(x,t;j, n) = (au*(x,f ; j,n), u*(x,< ;i,n)) (2.6) 

Let E 2 be divided into nonoverlapping rectangular regions (see Fig. 3(a)) referred 
to as conservation elements (CEs). As depicted in Figs. 3(c) and 3(d), the CE with 
its top-right (top-left) vertex being the mesh point (j, n) € fti is denoted by CE_(j, n) 
(CE+(i, n)). Obviously the boundary of CE_(j, n) (CE+(y, n)) is formed by subsets of 
SE(j, n) and SE(j - 1, n - 1) (SE(y + 1 , n - 1)). The current approximation of Eq. (2.1) 
is . 

f±(j,n) d =* h*-ds = 0 (2.7) 

Js(ce ;±<»)) 

for all (i,n) € Jli. In other words, the total flux leaving the boundary of any CE is 
zero. Note that the flux at any interface separating two neighboring CEs is calculated 
using the information from a single SE. As an example, the interface AC depicted in 
Figs. 3(c) and 3(d) is a subset of SE(j, n). Thus the flux at this interface is calculated 
using the information associated with SE(j,n). 

Because (i) The CEs associated with fij can fill any space-time region, and (ii) the 
surface integration across any interface separating two neighboring CEs is evaluated 
using the information from a single SE, the local conservation condition Eq. (2.7) leads 
to a global conservation relation, i.e., the total flux having the boundary of any space- 
time region that is the union of any combination of CEs will also vanish. 

With the aid of Eqs. (2.5)-(2.7), it can be shown that (see [1,8]) 

E± (j, ti)/ax = ±(1 - v 2 ) [(ti+ )" + (u+ )™g { ] + (1 T v) (vj - «"±x ) (2.8) 


where v — a At /ax is the Courant number, and (u+)" = (ax/2)(u*)". Note that here 
ax and At, respectively, represent the same mesh interval and time-step size which were 
denoted by ax/2 and At/ 2 in [1,8]. Using Eqs. (2.7) and (2.8), u* and (u+)?, which 
are considered as independent unknowns at the mesh point (j, n), can be solved for in 
terms of v"gl and («J)"±i if 1 — v 1 ^ 0. It is shown in [1,8] that, for all (j, n) € fli , 


q(j, n) = Q + q(j - 1, n - 1) + Q. q(j + 1, n - 1) (1 - 0) 

Here (i) q(j, n) is the column matrix formed by u? and («+)", and (ii) 

±v ±(1 — i/ 2 )' 

Tl 


(2.9) 


Q± = ( 1 / 2 ) 


/l±v ±(1 — 1 ^)\ 
\ Tl -l±») 


( 2 . 10 ) 


Eq. (2.9) defines a marching scheme. Because this scheme models Eq. (1.1) which is 
characterized by the parameter o, hereafter it will be referred to as the a scheme. 
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The a scheme is the only two-level explicit solver of Eq. (1.1) known to the authors 
to be neutrally stable, i.e., free from numerical dissipation. It also has the simplest 
stencil, i.e., a triangle with a vertex at the given time level and the other two vertices 
at the previous time level. Because the flux at an interface separating two neighboring 
CEs is evaluated using information of a single SE, no interpolation or extrapolation 
is required. Moreover, the a scheme is a two-way marching scheme, i.e., a backward 
marching scheme in which q(j,n) is determined in terms of q(j — 1, n + 1) and q(j + 
1, n+ 1) can also be derived from Eq. (2.8), the same relation that gives us the forward 
marching scheme Eq. (2.9). These and other nontraditional features of the a scheme 
are discussed in depth in [1,8]. 

In the above construction of the a scheme, we use the SEs and CEs of the mesh 
points marked by dots in Fig. 3(a). A similar construction can be performed by using 
the mesh points marked by crosses in Fig. 3(a). Let SI? denote the set of mesh points 
(j,n) in E? (crosses in Fig. 3(a)) where n = 0,±1,±2,±3,..., and, for each n, j = 
n, n ± 2, n ± 4, . . .. Let the SEs and CEs of Q 3 be defined by using Figs 3(b)-(d) with 
dots replaced by crosses. Obviously (i) the CEs of fl 3 also fill any space-time region, 
and (ii) the a scheme can also be constructed using the SEs and CEs of 0]. This new 
scheme is defined by Eq. (2.9) with (J, n) 6 £1%. 

Before we proceed further, some intricate points related to the above constructions 
will be clarified with the following remarks: 

(a) SE (i,n) may intersect SE(j',n') if (j, n) € fli and (j 1 , n') € Or- 

(b) Let (j,n) € 0i« Then (i) (j + l,n) € fi 3 , and (ii) CE+(j,n) and CE _ (j + l,n) 
represent the same rectangle in E?. However, because the function h* used in the 
evaluation of F+ (j, n) is tied to a pair of SEs associated with (2i, while that used 
in the evaluation of F-(j + l,n) is tied to another pair of SEs associated with 
flj, F+(j,n) = 0 and F_(j + l,n) = 0 represent two completely independent flux 
conservation conditions. 

To prepare for the development of the implicit solver to be described in the next 
section, we shall combine the above two independent schemes into a single scheme. 
The new scheme, referred to as the a(2) scheme, is defined by Eq. (2.9) with (j, n) € & 
where 12 1=F U fl 3 . Obviously, a solution of the a(2) scheme is formed by two 
decoupled solutions with each being associated with a mesh that is staggered in time. 
Several classical schemes also have this property. Among them are the Leapfrog, the 
DuFort-Frankel, and the Lax schemes [13]. 

We conclude this section with a review of other 1-D extensions [1,8] of the a 
scheme: 

(a) The a-p scheme is an explicit solver for Eq. (1.2) (/1 > 0). It reduces to the a 
scheme if p = 0. By using exactly the same procedure by which the a(2) scheme 
is formed from two independent a schemes, the a-p( 2) scheme is defined using two 
independent a-p schemes. 

(b) A scheme that has no numerical dissipation, such as the a scheme, generally can 
not be extended to solve the Euler equations. Hence, the a scheme is modified to 
become the a-c scheme. Stability of this scheme is limited by the CFL condition 
and 0 < € < 1 where t is a special parameter that controls numerical dissipation. 
If e = 0, the a-c scheme reduces to the a scheme. By using exactly the same 
procedure by which the a(2) scheme is formed from two independent a schemes, 
the a-e(2) scheme is defined using two independent a-c schemes. 

3. The implicit scheme 

An implicit solver for Eq. (1.2), referred to as the a-p(Il) scheme, will be dis- 
cussed in this section. Here “J* stands for “implicit”, and “1” is the identification 
number. This solver is the model for implicit time-dependent Navier-Stokes solvers 
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under development. It is constructed to meet two requirements given in the following 
discussion: 

(a) With a few exceptions, numerical dissipation generally appears in a numerical solu- 
tion of a time-marching problem. In other words, the numerical solution dissipates 
faster than the corresponding physical solution. For a nearly inviscid problem, 
e.g., flow at a large Reynolds number, this could be a serious difficulty because 
numerical dissipation may overwhelm physical dissipation and cause a complete 
distortion of solutions, lb avoid such a difficulty, the model solver is required 
to have the property that the numerical dissipation shall approach zero as the 
physical dissipation approaches zero. 

(b) The convection term and the diffusion term in Eq. (1.2) involve the spatial deriva- 
tives of first order and second order, respectively. Thus, in a spatial region where 
a solution is very smooth, the diffusion term is negligible compared with the con- 
vection term. As a result, the effective physical domain of dependence is more or 
less dictated by Eq. (1.1). To prevent excessive contamination of the solution by 
extraneous information, the implicit solver shall be required to become an explicit 
solver in the limiting case in which the diffusion term vanishes. 

Because of the requirements set forth in (a) and (b), the implicit solver will be 
constructed such that it reduces to the a(2) scheme if p = 0. The former differs from 
the latter only in the extra modeling involving the diffusion-related terms. Note that 
the presence of viscosity is felt through (i) the diffusion term — pd 2 u/dx 2 in Eq. (1.2), 
and (ii) the spatial diffusion flux component —pdufdx in the flu* vector 

h = (au — fidu/dz, u) (3.1) 

Note that Eq. (1.2) is the differential form of Eq. (2.1) if A is defined according to 
Eq. (3.1). Also, in this paper, any term in Eq. (1.2) or on the right side of Eq. (3.1) is 
considered to be convection-related if it is not diffusion-related. 

t,n 


j = 0 j = J 



Figure 4. — The space-time mesh for the implict solvers (J = 6). 

To construct the a-ft(Il) scheme, consider the mesh depicted in Fig. 4 (J > 4). 
We assume that (i) u = uj(a:) at t = 0, (ii) u = tit(f) at x = 0, and (iii) u — un(t ) at 
x = Jax, where uj(ar), «£,(<), and ur ( t) are given functions. Moreover, for the current 
case, (i) Qi and &2 are restricted by the conditions n > 0 and J > j > 0, (ii) CE±(j, n) 
are not defined if n = 0, (iii) CE_(j, n) is not defined if j = 0, and (iv) CE + (j,n) is 
not defined if j = J. Items (iii) and (iv) imply that only one conservation condition 
is associated with a boundary mesh point. Obviously, the definition of SE(j, n) also 
needs to be appropriately modified if j = 0, or j = J, or n = 0. 

Eq. (2.2) will still be assumed. We also assume that, for n = 0, 1, 2, . . ., 

«o =«£(*") (««)o «3 = u*(f n ) (u t )S = «*(*") (3-2) 
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where u^(i) dui(t)fdt and ur(1) = dun(t)/dt. Thus, only one unknown, i.e., (ti*)" , 
and one conservation condition are associated with a boundary mesh point (j, n). 
Furthermore, for an interior mesh point (j, n), we replace Eq. (2.4) with 

(«.)? = -«(«,)? + £ [(«,)J +1 - («,)?-,] J>>>0 (3.3) 

Eq. (3.3) is the numerical analogue of the differential condition Eq. (1.2). A comparison 
between Eqs. (1.2) and (3.3) reveals that the diffusion term —pfPu/dx 2 at an interior 
mesh point € fli (H 2 ) is modeled by a central-difference approximation involving the 
values of u r at two neighboring mesh points £ Q 2 (fti). Eqs. (2.2) and (3.3) imply 
that, for J > j > 0 and (*,f) 6 SEO', n), 

u*(i,t;i,n) = «" + («.)"(* -*j) + { 2 “ [( u *)?+i ~ K)"-i] -“(«*)? }(*"*") ( 3 - 4 ) 
Next, as a result of Eq. (3.1), Eq. (2.6) is replaced by 


h*(x,t;j,n) = (au*(x,t;j,n) — /x«*(x,t ;y,n), u*(x,t ;j,n)) 
Here, for any (x,t) € SE(j, n), 


„„ f [(«*)" + («*)" +1 ] A ift>t n ; 

l ((«.);+(«*)"" ] a **<*"• 


(3.5) 


(3.6) 


Consider any point (x, t) on the line segment joining the mesh points (j, n) and 
(j,n— 1). Then (x,t) belongs to both SE(j, n) and SE(j,n— 1). According to Eq. (3.6), 
for this point (x,f), 

«£(*,< *, j,n) = ul(x,t ]j,n - 1) = [(«*)" + (t»*)" -1 ] /2 (3.7) 

Thus, the same numerical diffusion flux component is assigned to the point (x,f) re- 
gardless of whether it is considered as a point in SE(J, n) or a point in SE(j , n — 1). 

With the above modifications, the a-/i( 71) scheme is defined by assuming Eq. (2.7). 
Note that: 

(a) At a mesh point € fli(^2)> the diffusion-related terms in E<p. (1.2) and (3.1) are 
modeled using interpolations that may involve the numerical values of the mesh 
points € fts(Qi). This contrasts sharply with the modeling of the convection- 
related terms which uses no interpolation. 

(b) In the a(2) scheme, the two sets of numerical variables associated with fii and Sli 
are completely decoupled from each other. Contrarily, they are “glued” together 
in the a- ft (II) scheme through the interpolations referred to in (a). 

To proceed, let a = ft At /(ax) 2 , and (ti*)" 1= (At/2)(ti«)". Also let 


(S + )J « (1 - ..)«?„ + (5 - 1W«J)? - (1 - - «)(<•!)"* 1 


- TrOi )“« (n = 0, 1. 2, ... J = 0, 1, 2, .... J - 2) 


(3.8) 


(s+»-i = (1 - *-)«” - «(«i)3-. - (i - <*)(«:» - K«, + )5 

(S-); « (1 + +(j+ i)a(«i)? + 0 

- y («J )"-J (n = 0, 1, S , . . . ; j = 2, 3, 4, . . 


(n = 0, 1,2, ...) (3.9) 


)(<)?-! 


‘ J 


J)) 


(3.10) 
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(S'-)" = (1 + »/)«o + <*(«*)" + (1 — <*)(“* )o + K u ?)o (n = 0,1,2,...) (3.11) 

By using Eqs. (2.2), (3.2), (3.4)-(3.6), and (3.8)-(3.11), Eq. (2.7) implies that 


(1 + v)u*j ± (1 - v 2 + o)(w+)7 + (^ =F l)a(«i")"±i - y(«J)” t i 
= (5*);- 1 (n = 1, 2, 3, . . . ; j = 1, 2, . . . , J - 1) 

(1 + a)(«tf£ - «(«+)" = (5 + )r 1 - (1 - u)uZ - *(«+)S (n = 1,2, . . .) 

<*£ £-1 - (1 + «)(«?)" = (5-)r 1 - (1 + + K«* + )3 (n = 1,2, .. .) 

Eqs. (3.12)-(3.14) define the a-p(/l) scheme. 

Note that Eq. (3.12) represents a pair of equations for each (j, n). Let 1 — i^jtO. 
Then these equations are equivalent to 

-cr(«+)7_ 1 +2(l-^+a)(u+);-a(«+); +1 = (1 +»)(S + )?- 1 -(l-v)(S-)?- 1 (3.15) 
and 

«? = *{( 5 +)r 1 + + « [(■*»+» - («? )"-ii } (3i6) 

where j = 1,2,3,...,/ — 1. In the following discussion, Eqs. (3.12)— (3.14) will be 
replaced by Eqs. (3.13) — (3.16). 

Let the marching variables at the (n — l)th time level be given. With the aid of 
Eq. (3.2), the expressions on the right sides of Eqs. (3.13)-(3.15) can be considered as 
given source terms. Thus these equations form a tridiagonal system of J + 1 equations 
for the 7+1 unknowns (uj )?, j s 0, 1, 2, . . . , J. It can be shown that [3] the coefficient 
matrix is strictly diagonally dominant in rows and columns, i.e., stability of the Thomas 
algorithm [13, p.99] is assured, if 

v 2 < 1 and a > 0 (3.17) 


(3.12) 

(3.13) 

(3.14) 


Upon obtaining (u+ )” , j = 0, 1, 2, . . . , J, the other unknowns ti", j = 1, 2, . . . , J — 1, 
can be obtained using Eq. (3.16). 

Using a procedure described in [1,8], the stability of the a-p(il) scheme was studied 
using the von Neumann analysis [3]. Because there are two independent marching 
variables at each (j, n), the a-/i(/l) scheme has two amplification factors. The principal 
amplification factor A + and the spurious amplification factor A_ are given by 


A± = 


-Q ± y^ 2 + (1 - »*)* - q 2 (l - cos#) 2 
1 — V 3 + o(l — cos 0 ) 


where 0, —v < 0 < v, is the phase angle variation per ax, and 


Q = f a( 1 — cos 6 — j/ J sin 2 0) + »V(1 — u 2 ) sin# (» 1= >/— T) 


(3.18) 


(3.19) 


It is concluded from numerical evaluations of Eq. (3.18) that |A±| < 1, i.e., the a-/i(/l) 
scheme is stable, if Eq. (3.17) is satisfied. 

Note that many other implicit solvers are unconditionally stable. However, the 
price paid for this “desirable” property usually is excessive numerical dissipation. More- 
over, the use of a time-step size that is greater than that allowed by Eq. (3.17) generally 
results in a less accurate time-dependent solution. Thus we do not consider the more 
restrictive stability condition Eq. (3.17) to be a disadvantage of the a-p(Jl) scheme. 
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One can also conclude from Eq. (3.18) that: (i) A± are reduced to the amplification 
factors of the Leapfrog scheme if p = 0 and (ii) A + is reduced to the amplification factor 
of the Crank-Nicolson scheme [13, P-112] and A. = — 1 if a = 0. 

The fact that the amplification factors of the a-p( 71) scheme are related to those 
of two celebrated classical schemes is only one among a string of similar coincidences . 
Other coincidences are summarized in the following remarks [1,3,8]: 

(a) Because the a-p( 71) scheme reduces to the a(2) scheme if /i = 0, the amplification 
factors of the a(2) scheme are those of the Leapfrog scheme . 

(b) The amplification factors of the a-p(2) scheme reduce to those of the Leapfrog 
scheme if fi = 0, and to those of DuFort-Frankel scheme if a = 0. 

(c) If € = 0, the amplification factors of the a-c(2) scheme reduce to those of the 
Leapfrog scheme . On the other hand , if e = 1, these two factors become the same 
function of the Courant number and the phase angle . It so happens that this 
function is also the amplification factor of the highly diffusive Lax scheme. 

In summary, the amplification factors of the Leapfrog, the Crank-Nicolson , the DuFort- 
Frankel, and the Lax schemes are related to those of the current model schemes. 

In [3] , by using the von Neumann analysis, it is shown that the a-p(71) scheme is 
more accurate if a spurious component is filtered out from its initial conditions. Let 
it/(x) = dui(z)/dz where u/(x) is the given function introduced earlier. Then it is 
shown in [3] that, after filtering, the initial conditions are (i) 

ttj = {«/(* i+ i) + «/(*,-_ i) + 2 u/(xy) + (ax/ 2) [«/(*,_!) - *j(*,+i)]}/4 (3.20) 

(«+)“ = («/(x i+ i) - + (ax/ 2) [2ti/(x ; ) - u/(x i+ i) - i/fo-!)]}/* (3.21) 

where j = 1, 2, 3, .... J - 1, (ii) tig = u/(x 0 ) and ug = tij(xj), and (iii) 

(«*)o = {«/(*i)“ «/(x 0 ) + (ax/2) («/(xo) - «/(*i)]}/2 (3.22) 

(t tt)°j = {«/(x,) - ii/(*.r_i) + (ax/2) [uj(xj) - ti/(x/_i)]} /2 (3.23) 

This section is concluded with a brief discussion of another implicit scheme, re- 
ferred to as the a-/i(/2) scheme. It differs from the a-/i(Il) scheme only in one aspect, 
i.e., Eq. (3.3) is replaced by 

(«<)" = -«(«•)" + + *?-i - *9) J >i>° (3-24) 

As a result, Eq. (3.4) must be modified accordingly. The details of this scheme will be 
given in [3]. Note that (i) stability of the a-/i(/2) scheme is also limited by Eq. (3.17), 
and (ii) the a-/i(/l) and a-/i(J2) schemes are identical if p = 0 or a = 0. 

4. Numerical Results 

Three test problems will be used to evaluate the accuracy of the a-fi(1 1) scheme. 
In the first problem, we consider a special case of Eq. (1.2) (a = 0 and ft = 1) with 
0 < x < 1 and t > 0. The analytical solution u = U(x,t) with the initial/boundary 
conditions (i) u = 0, x = 0 and 1, t > 0; (ii) u = 2x, 0 < x < 0.5, t = 0; and (iii) 
ii = 2(1 — x), 0.5 < x < 1, t = 0, is given on p.15 of [14]. Obviously, U( 0.5 + x',t) = 
Z/(0.5 — x',<) for lx*! < 0.5. 

The numerical solutions given in Table 1 are generated assuming ax = 0.1, and 
At = 0.01. Note that the a-p(Jl) scheme, like the Crank-Nicolson scheme, is uncondi- 
tionally stable in this case (a = 0). From this table, one concludes that the numerical 
solutions from these two schemes become closer as time increases. This conclusion is 
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consistent with the fact that the principal amplification factor of the a-/i(/l) scheme 
is identical to the amplification factor of the Crank-Nicolson scheme. 


TABLE 1 .—COMPARISON OF THE SOLUTIONS OF THE FIRST TEST 
PROBLEM OBTAINED USING THE ANALYTICAL METHOD, THE 
a-|x(ll) SCHEME AND THE CRANK-NICOLSON SCHEME 



II 

X 

0.1 

0.2 

0.3 

0.4 

wm 


analytical 

0.1996 

0.3966 


0.7201 

0.7743 

t = 0.01 

a-fi(II) 

0.1992 



0.7286 

0.7768 


Crank-Nicolson 

0.1989 



0.7381 



analytical 

0.0934 


0.2444 

0.2873 

0.3021 

t = 0.1 

a-fKIl) 

0.0947 

BEa 

0.2481 

0.2917 



Crank-Nicolson 

0.0948 

BE! 

0.2482 

0.2918 



In the sceond problem, we consider a special case of Eq. (1.2) (a = ft = 1) with 
0 < x < 1 and f > 0. The steady-state solution with the boundary conditions (i) u = 1, 
* = 0; and (ii) u = 0, x = 1, is given on p.155 of [13]. 

Let (i) «/(*) = 1 — x, (ii) ax = 0.25, and (iii) v = 0.5 (i.e., At = 0.125). After 40 
time steps, the numerical values of u at x = 0.25, 0.5, and 0.75 are 0.8356, 0.6234, and 
0.3504, respectively. The corresponding exact steady-state values are 0.8347, 0.6225, 
and 0.3499, respectively. The numerical results are accurate to almost three significant 
figures even though only three interior mesh points are used. 

In the third problem, again we consider a special case of Eq. (1.2) (a = 1 and 
p = 0.1) with 0 < x < 1 and t > 0. The functions ti/(x), tit(t) and ttjr(i) are defined 
such that they are consistent with a special solution to Eq. (1.2), i.e., 

ti = u«(x, t) d = exp(— 4 x 2 /it) sin [2x(x — at)] (4.1) 


Let u e *(x,i) = f du e (x,t)/dx. At any time t = t n , let 

J - 1 


* (/-DexpMxV") § |tt “ " tle(xj,<n)l 


(4.2) 


iiK) d = 


(J + 1)2 


xexp(-4rV n ) ^ ~ 


(4.3) 


L\(u) and Li(u*) are two error norms (per mesh point) which are normalized by the 
decay factors of u e (x,i) and u es (x,t), respectively. 

Let J = 80 (i.e., ax = 1/80) and v — 0.8 (i.e., a t = 0.01). Then Li(u) = 
0.2312 x 10 -2 and Li(tx*) = 0.1363 x 10 -1 at t = 1 (i.e., n = 100). Through numerical 
experiments, it has been shown that both Li(u) and Li(tt r ) at a given time t are 
reduced by a factor of 4 if both ax and At rue reduced by half. 

5. Conclusions and Discussions 


In the explicit a-/i scheme [1,8], the diffusion term in Eq. (1.2) is not modeled, 
i.e., Eq. (2.4) is assumed. Also the diffusion term in Eq. (3.1) is modeled with no 
interpolation or extrapolation. Obviously, it can be used only when the diffusion term 
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is small compared with the convection term. Contrarily, the diffusion terms in both 
Eqs. (1.2) and (3.1) are modeled in the current implicit solvers. 

A more complete discussion of the a-p(Il) and a-/i(/2) schemes will be given in 
[3J. In this paper, we shall (i) study the stability, dissipation, dispersion, consistency, 
and truncation error of the above two schemes, and (ii) discuss other explicit solvers 
for Eq. (1.2) in which the diffusion terms in both Eqs. (1.2) and (3.1) are modeled. 
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